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Abstract 

In this paper we revisit the problem of Brownian motion in a tilted periodic potential. We 
use homogenization theory to derive general formulas for the effective velocity and the effective 
diffusion tensor that are valid for arbitrary tilts. Furthermore, we obtain power series expansions 
for the velocity and the diffusion coefficient as functions of the external forcing. Thus, we provide 
systematic corrections to Einstein's formula and to linear response theory. Our theoretical results 
are supported by extensive numerical simulations. For our numerical experiments we use a novel 
spectral numerical method that leads to a very efficient and accurate calculation of the effective 
velocity and the effective diffusion tensor. 

1 Introduction 

Brownian motion in periodic potentials is one of standard models in condensed matter physics. Ap- 
plications include the modeling of Josephson junctions, polymer dynamics, superionic conduction, 
dielectric relaxation, plasma physics and surface diffusion [lj. A detailed discussion and extensive 
bibliography can be found in O [27] . 

The goal of this paper is to study Brownian motion in a tilted periodic potential for arbitrary 
values of the drift and of the tilt (external forcing). The dynamics of the Brownian particle is 
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governed by the Langevin equation 



q = -V Q V(q) + F- 1 q + v / 2 1 p~ 1 W, (1) 

where V is a periodic potential with period L (in each direction), F denotes a constant external 
force, so that the effective potential is 

V eS (q) = V(q)-F.q, (2) 

7 is the friction coefficient, j3 is the inverse temperature and W(t) is a standard Brownian motion 



on 



The main objective is the calculation of the drift and diffusion coefficients which are defined as 

v = lim w) -»(»)) (3) 

t—too t 

^ D = Hm (q(t) - (q(t))) ® (Q(t) - (q(t)))) _ (4) 

t—too 2i 

Here (•) denotes ensemble average and ® stands for the tensor product. Explicit formulas for these 
coefficients are available only in the overdamped limit and mostly in one dimension. An exact 
analytical formula for the effective velocity of an overdamped Brownian particle moving in a one 
dimensional tilted periodic potential was obtained many years ago by Stratonovich ([29 ,[30, Ch. 
9]). A corresponding analytical formula for the diffusion coefficient was obtained and analyzed 
more recently |25 |, 1241 [22], an d verified in an experimental realization of the model involving rotat- 
ing optical tweezers [9|. Simpler algebraic formulas were deduced from these for the special case 
of piecewise linear potentials in [13J. Only potentials with very specific geometries can lead to 
analytical formulas in dimension higher than one [H [3] . A wealth of information on the problem of 
Brownian motion in a tilted periodic potential in one dimension can be found in \27\ Ch. 11]. 

It is well known that the equilibrium diffusion coefficient (i.e., the diffusion coefficient in the 
absence of an external drift) and the drift, or, rather, the mobility are related through the famous 
Einstein formula: 

D = /rV, (5) 



with 



a = lim VfU 

\F\->0 



The validity of this formula has been proved rigorously for several models [19] , including that of a 
Brownian particle in a tilted periodic potential [2H]. Formulas of the form §5§ can be understood 
in the more general framework of linear response theory and of the Green-Kubo formalism |26l [18] . 
A recent rigorous analysis of the Green-Kubo formalism for the calculation of the shear viscosity 
can be found in [16] . 

The main goal of the present paper is to investigate the validity and usefulness of corrections 
to linear response theory. In particular, we calculate all terms in the power series expansions 
(with respect to the forcing F) for the drift and diffusion coefficients and we use these in order to 
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calculate corrections to Einstein's formula ([5]). Our analysis is based on the formalism of averaging 
and homogenization [23J. From this formalism we know that both drift and diffusion coefficients 



can be expressed in terms of the solution of appropriate Poisson equations (14)-(17). Details are 
presented in the next section. 

For simplicity of notation and presentation, we will restrict our calculations for corrections to 
linear response theory to the one dimensional case d = 1, and hence hereafter drop vector and 
tensor notation. Completely analogous formulas are applicable in multiple dimensions. We present 
our results in detail in Section [3j and summarize them here rather imprecisely: 

Proposition 1.1. The drift and diffusion coefficients admit the asymptotic expansions 

U = Y,F% , (6) 



£>1 



and 

D = 



i>0 



n=l 



(7a) 



The coefficients Vj, £ n ,-, n, j = 1,2,... can be computed in terms of solutions to Poisson 
equations for the generator of the equilibrium dynamics F = 0. In particular, the higher order 
corrections to the drift and diffusion coefficients are not compatible with an extension of the Einstein 
relation |5p beyond the linear response regime F — > 0. 

For the case of a symmetric potential (V{q) = V(—q)), then V n = for even n and T, n £ = 

= for odd £. 

Thus, it is possible, in principle, to calculate the drift and diffusion coefficients of the nonequi- 
librium dynamics ([T]) in terms of the equilibrium dynamics 

Q = -V'{q) - 75 + V^P^W (8) 

for at least a finite interval of values of F. 

The validity and usefulness of the power series expansions ^ is tested by performing numerical 
experiments. For the calculation of drift and diffusion coefficients we need to solve Poisson equations 
of the form 

-C4> = f(p,q), (9) 

with C being the generator of the Markov process {q, p} with p = q. We solve equations of the 
form (|9]| using a spectral method 12 U that is an extension of Risken's continued fraction expansion 
method [27]. By comparing the results obtained using our spectral method with results obtained 
from (the computationally more expensive) Monte Carlo simulations, we demonstrate that our 
method performs very well. 
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The rest of the paper is organized as follows. In Section [2] we present the formulas for the 
drift and diffusion coefficients obtained using homogenization theory. In Section [3] we calculate the 
power series expansions for the drift and the diffusion coefficient. In Section [4] we present results of 
numerical simulations on the calculation of U and D. Section [5] summarizes our conclusions. The 
details of the spectral method for the solution to the Poisson equation are presented in Appendix [A} 
Some discussion of how our formulas relate to an alternative approach developed in [2] can be found 
in Appendix [B] 

2 Formulas for the Drift and Diffusion Coefficients 

We start by writing as a first order system, in d = 1 dimension: 

q = p, p = -V'(q) + F-jp+ sfinP^W. (10) 
The process {q, p} is a Markov process with generator 

c = P -d g + {-d q v + f) ■ d p + 7 ( - P • d p + p~ l dl) . (ii) 

The Fokker-Planck operator, i.e. the L 2 -adjoint of the generator, is 

C* = -p ■ d g + (8 g V — F) ■ d p + jd p ■ [p + /r%) . (12) 

The potential function V has period L. We can use homogenization theory [22\ [T2l [20] to prove 
that the rescaled process 

<f(t):=eq(t/e 2 )-— , (13) 
e 

where U is the effective drift as defined below, converges weakly on C([0,T];M) to a Brownian 
motion with diffusion coefficient D. To write down the formulas for the drift U and the diffusion 
coefficient D we need to consider the process X(t) = {q(t), p(t)) defined on X := T x 1 where T 
denotes a one-dimensional circle with length L corresponding to the period of the potential V. The 



generator and Fokker-Planck operator of this process are still given by formulas (11) and (12) but 
now restricted on X, with periodic boundary conditions with respect to q. It can be shown |28j 
that X(t) is an ergodic Markov process with invariant measure ^^{dpdq) that has a smooth density 
PpiPiQ) with respect to the Lebesgue measure on X. The invariant density is the unique solution 
of the stationary Fokker-Planck equation on X: 

C* PP = 0. (14) 

The drift is then given by the average of the momentum with respect to pp over X: 

U = / PP/3(p,q)dpdq. (15) 

The diffusion coefficient is given by the formula 

D = / (p -U)(j)pp{p,q)dpdq (16a) 

= / ' {d P <t>? pp(p,q)dpdq, (16b) 
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where <j) is the solution of the Poisson equation 



Ccj) = p-U, 



dpdq = 0. 



(17) 



Equations (14) and (17) are equipped with periodic boundary conditions in q and suitable inte- 



grability conditions \22 \ [12 | [20] . Formula (16b), which shows that the effective diffusion tensor is 



positive semidefinite, follows from (16a) after an integration by parts. 



It is possible to prove that both U and D are analytic functions of the forcing F. This has 
been proved for the drift in [6j (in fact, in this paper the analyticity of the drift with respect to 
the forcing is proved for several models including systems of coupled Fokker-Planck equations). A 
similar analysis can be used to prove the analytic dependence of D on F. 



3 Corrections to Linear Response Theory 

In this section we solve perturbatively equations ( Jl~4| ) and ( |17| ) in one dimension, in order to obtain 
the power series expansions (|6]). Calculations of this form are quite standard when investigating 
the effect of colored noise on the drift and diffusion coefficients, e.g. [T3JIB1E2]. Related calculations 
have presented recently in |16j . 

The main result of this section is a precise formulation of Proposition |1.1| To state the re- 
sult, we need to introduce some notation. We denote by Hq the Hamiltonian of the unperturbed 
(equilibrium) dynamics ([!]): 

H (p,q) = \p 2 + V{q). 
The invariant density of the unperturbed dynamics on X is denoted by p: 

p{q,p) = \e-^ Ho ^ p \ Z = [ e-^ Ha ^ dqdp. (18) 
% Jx 

We will work in the weighted 1? space H = L 2 (T xR;p), The inner product in this Hilbert space 
will be denoted by (•,•)/?• The generator of the unperturbed dynamics can be written in the form 

£ = A + jS, (19) 

where 

A = pdq — d q Vd p and S = — pd p + /3 -1 <9p, 

denote the reversible and irreversible parts, respectively. The operators A and S are antisymmetric 
and symmetric, respectively, in 7-L. We introduce now the creation and annihilation operators |31| 

EZ1 CEDE! 

a + := —d p + /3p and a" := d p . (20) 
These two operators are %-adjoint: 

{a + f,h)p = {f,a-h)p, Vf,hen. 
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Proposition 3.1. The drift and diffusion coefficients admit the asymptotic expansions 



(21) 



e>i 



and 



D 



£>0 



ni 



n=l 



p~'% + Y. f 'Y. : 



£>1 n =l 

The coefficients Vn, E^, S^t k < £ = 1,2, . . . are given by the formulas 

Vi = j feppdpdq = (3 J fa^xppdpdq, 

= / P4>l-nfnPdpdq 

^ni = fi~ X \ 4>i- n dpfnpdpdq 
where /,-, (ftj, j = 0, . . . are solutions to the (adjoint) Poisson equations 

-C fj = a + /j-i, /o = 1, fjp = 0, j = l,2,. 



4> p = o, 



-c 



09 j 



» / <l>jP = ~ E / /r^J- 



-P 3 



(22a) 
(22b) 

(23a) 
(23b) 
(23c) 

(24a) 
(24b) 

(24c) 



We have used the notation Co = —A + 75 to denote the %-adjoint of the generator £q. 



Remark 3.1. The expansion formulas for the drift and velocity are consistent with the exact 
statistical reflection symmetry q — > —q, p — > —p, and F — > —F in the stochastic system (10) or 

V(—q). Since the drift is 
when i is 



infinitesmal generator (11) when the potential is symmetric: V(q) 
odd and the diffusivity even under reflection, this implies that the coefficients Vt 
even and the coefficients Ti n i = H n £ = when I is odd. One can verify that our formulas do 
indeed have these vanishing properties under symmetry of the potential, noting that the operators 
a + and a~ are odd under reflection, whereas £$ and Cq are even. By uniqueness of solutions of the 



Poisson equations ( 24 ) and the symmetry properties of the operators, inhomogeneity, and auxiliary 



conditions, we first verify by induction that the functions fj have even reflection symmetry for even 
j and odd reflection symmetry for odd j. Then we similarly induce that the functions 4>j have 
odd reflection symmetry for even j and even reflection symmetry for odd j. Finally, p (18) has 



manifestly even symmetry under reflection symmetry. Therefore, when t is even, Vi can be checked 
to be the periodic integral of an odd function and when I is odd, E n £ and are periodic integrals 
of odd functions, and so vanish. 
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Using the notation that we have introduced in this section, Einstein's formula (linear response 
theory) can be written in the form 



However, formula (22a) shows that it is not true that a similar simple relation holds for higher 
order terms in the expansions for the drift and the diffusion coefficients. In particular, it is not 
true that 

D n = P~ 1 V n+ i, n = l,..., (25) 



but instead there is a non-trivial correction to (25) that is given by the second term on the right 



hand side of (22a). As an example, we present the formula for the diffusion coefficient that is valid 
up to 0(F 3 ): 

D = (3- 1 V 1 + F^- 1 V 2 + Jpfofipdpdq} (26) 

+ F 2 ((3- 1 V 3 + [ p<p 1 f 1 pdpdq + [ P M2pdpdq)+0(F 3 ). 

v Jx Jx ' 

Notice that the calculation of the next two terms in the expansion for the diffusion coefficient 
requires the solution of an additional Poisson equation, in order to compute </>i, as well as the 
calculation of three additional integrals. 

Similarly, it is not true that the Einstein relation ^ can be extended away from F = in 
the form: 



because of the presence of correction terms in Eq. (22b). This issue is investigated numerically in 
the next section, see Figures |4| and [5} The relation Eq. ([27} was indeed hypothesized in [7] , but |24j 
showed through analytical and numerical studies that while it seems qualitatively correct, and is 
quantitatively correct in the three limits F — > 0, F — > oo, and (3 — > oo, it is not quantitatively 



accurate for general parameter values. Our results in Proposition 3.1 give quantitative formulas for 
this discrepancy, for example, through third order: 

D = r ld ^p- + Fr X J fodpfipdpdq (28) 

+ F 2 $- l { [ c/ )0 d p f 2 pdpdq+ [ fadpftpdpdq). (29) 

K Jx Jx J 

The violation of the Einstein relation for F ^ in the model under consideration, and other 
nonequilibrium steady-state models, was recently analyzed by |2J from a different nonperturbative 
perspective, expressing the correction terms with respect to various time-correlation functions of 
the dynamics. But as we discuss in Appendix [B] our framework based on perturbation expansions 
of the equations from homogenization theory appear to yield more easily computable expressions. 
We remark also that [T7] have examined deviations from the Einstein relation in the context of 
stochastic tracer dynamics in a random environment. 
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Proof of Prop. 3.1 We start with the analysis of the stationary Fokker-Planck equation (14). 
We set 

pp(p, q) = Piv, q)f(p, q) (30) 



We substitute (30) into (14) and use the symmetry and antisymmetry of S and A, respectively as 
well as equation (20) to obtain 



Co + Fa+)f = 



(31) 



Equation (31) is posed on X := T x R and is equipped with periodic boundary conditions with 
respect to q as well as the normalization condition 

fpdpdq = 1. 



We look for a solution to (31) as a power series expansion in F: 



N 



f(p,q) = Y J FJ fj(Pi 

3=0 



(32) 



The normalization condition becomes 

N 



y2 FJ //(Pi q)p(p> q) d P d q = L 



3=0 

This condition has to be satisfied for all F E M. which implies that the following normalization 
conditions should be satisfied 



{f Q ,l)p = l, (/,-, 1)^ = 0, J = 1,2, . . . 



We substitute the expansion for / into (31) to obtain the sequence of equation 



£o/o 

-£ofj 



0, 

a + fj-i, j = 1,2, 



(33) 



(34a) 
(34b) 



The above equations are of the form 



— £oV' = u - (35) 

The null space of Co, as well as its %-adjoint Cq is one-dimensional and consists of constants. 
Consequently, the solvability condition for equations of the form ( 35 ) is that 

(M)/3 = o. 



(36) 



Provided that the solvability condition (36) is satisfied, the Poisson equation (35) has a unique 
mean zero solution, (ip, 1)^ = 0. We correspondingly define the operator (— £o) _1 on the subspace 
of functions satisfying (36) to be this unique mean zero solution. 



From the first equation in ( 34 ) and the normalization condition we deduce that 

/o = l- 



(37) 
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The properties of the operators immediately yield that the solvability condition is satisfied for 
all equations for fj , j = 1 , . . . : 

(a + fj-x,l)i3 = (fj-i,a~l)p = 0. 



The solution of Equations ( 34b ) can be written as 

fj = {-£o)~ 1 a + fj-i = fcfj-i, 
where K, is the %-adjoint of K, := a~ (— Cq 1 ). Consequently, 

fi=tn, j = 0,1,... (38) 
Thus, we have obtained a power series expansion for the invariant distribution in powers of F: 



p(p,q;F)=p(p,q) 1 + ^F^'l 
from which we immediately deduce the expansion for the effective drift: 

v = J2 pj (P'fjh 

i>i 
i>i 

i>i 

= r^F^i^-zcn)^. 

i>i 



(39) 



In particular: 



Vj = /3 -1 (l,a~/C J l)^. 



(40) 



(41) 



Now we proceed with the analysis of the Poisson equation (17) which, in view of (40), rt30b, (32) 



and ( 37 ) , can be written as 



(42) 



with Vj given by (41). The generator of the perturbed dynamics is 



£ = Co + Fa 



where Cq is given by (19). We look for a solution of (42) in the form of a power series expansion 
in F: 

4>{p,q) = X^WP'"?)- 

j>0 
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We substitute this expansion into Equation (42 ) to obtain the sequence of equations (recalling from 
Eq. d37b that f = 1): 



£o4>o = P, (00, = 



r=l 



, <Pj-r)p J 



1,2, 



(43a) 
(43b) 



Equation (43a) is precisely the Poisson equation for the unperturbed dynamics F = 0. Now we 
show that the solvability condition (36) is satisfied for equations (43b). We need to show that 

^• = (a~^-i,l>/3- (44) 



Lemma 3.1. The solvability condition (44) is satisfied for all j > 1, and moreover the relation 

(a - </>o, fk)p = {a~<Pm, fk-m)p fork = m,m + l,... (45) 

holds for all m > 0. 



Proof. Our strategy pivots on the observation that if we can establish (45) for m = 0, 1,2, . . . , M, 
then the solvability condition (44) follows for j = l,2,...,M + l: 

= /?" 1 (a~(-£ )~ 1 « + l,^ J_1 l>/3 = (a~<t>o,fj-x)p = (a~<pj-i, / )^ 
= {a~<j>j-i, l>/3, 



using Eq. ( 45 ) with m = j — 1 in the penultimate equality. 

We now proceed to establish (45) inductively for m = 0, 1,2, . . .. The case m = is trivial. 
Suppose now Eq. (45) has been shown for all m = 0, 1,2, ... , M; we will show that (45) also follows 
for m = M + 1. To this end, it is useful to introduce the operator P which projects orthogonally 
onto the hyperplane in % orthogonal to constants: 

Pg:=g-(l,g) p . 



Since, by the above argument and the induction hypothesis, the solvability condition for (43b) is 
satisfied for j = 1, 2, . . . , M + 1, we can write 

M+l 

4>M+1 = {-C-oT l Par (j) M - ^2 {fr,4>M+l-r)p, 



r=l 



where the second sum of constants is included to meet the side condition in Eq. (17), as we have 
defined (— £o) -1 to yield a mean zero solution. But the operator a~ will kill these constants, and 
therefore, for k = M + 1, M + 2, . . ., we can write: 

(a~(f)M+i,fk-M~i)p = (a~ (-Co^fa^ 4> M , fk-M-i)p 

= (fCPa~<pM,fk-M-i)p = (a~(/>M,P^fk-M-l)p 

= (o _ 0A/,A^/fe-Af-l)/3 = (a~(f>M, fk-M)p 

= {a~<t>o,fk)p- 
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In the penultimate equality, we used PAC = K. from the fact that by definition of (— Cq) , 



-Cr 



Cq) ; the final equality follows from the induction hypothesis. 



□ 



Now we derive (22a). Using the centering condition in (17) we have that the diffusion coefficient 
is given by 



D 



p4> p/3 dpdq 

J2J2 Fn+e [Pfafnpdpdq 
£>0n>0 ^ X 
oo r „ 

^2 ^2 pr \ p r - n dpdq 

r=0 n=0 ^ X 

oo T „ r 

V] F r p 4> r f p dpdq + S2 / P <p r ~n fnP dpdq 

r=0 I/* n=l J * 

oo j~ r „ 

/3~ 1 (cj) r ,a + l)i3 + '^2 / P <t> r -n fnP dpdq 

n=l ^ X 

r 



r=0 



r=0 



n=l 



yielding (22a 



We can also alternatively restructure this expansion as follows, using the relations (44) and 



(45): 



oo r ,. 

D = ^ ^2 ^ \ P ^ r - n dpd<1 
r=0 n=0 ^ X 

oo r 

r=0 n=0 

oo r 

= P' 1 ^2 ^ ^2 [( a_ ^r-n., fn)/3 + (<Pr-n, a~ f n ) ] 
r=0 n=0 

oo r 

= /T 1 J2 J2 i^'^' MP + &r-n, arf n )p] 
r=0 n=0 
oo r 

= p- 1 Y. FT 12\y^^^-n,a-u)p\ 



r=0 n=0 

oo 



r=0 

d£7 



+ 1)1^+1 + E(0 r _„, a / n ) 



ra=0 



r=l n=l 



establishing the statement ( 22b ) in the proposition. In the last equality, we used that a fo = 0. □ 
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0.2 



0.4 



0.6 



F/F 



0.8 



1.2 



Figure 1: Solid Lines: U as a function of F for 7 = 0.01, 0.1 and 1. Markers: Monte Carlo 
estimates. Other parameters of the simulations are Vo = 7r 2 /16, /3 = 1.2V^~ , and L = 2ir. Following 
convention, the force variable is scaled by the critical force F c ~ 3.367\/Vo at which the effective 
potential Q becomes monotonic, and the drift is scaled by the value Ul = F/j it would have in 
absence of the periodic potential <f>(x). For the Monte Carlo simulations we use an Euler-Maruyama 
scheme with a time step At = 0.1 integrating over N = 5000000 time steps (after 100000 time steps 
of a transient integration interval) and averaging over 5000 trajectories. 

4 Numerical Simulations 



In this section we present results of numerical simulations that corroborate the theoretical results 
presented in the previous section. The calculation of the drift and diffusion coefficients is based 



on the numerical solution of the hypoelliptic boundary value problems (14) and (17) as well as the 



calculation of the integrals (15) and (16a). Both PDEs are solved using a spectral method that 



relies on the expansion of the solution of the stationary Fokker-Planck and the Poisson equations in 
a Fourier-Hermite expansion. This method is adapted from Risken's continued fraction expansion 
method [27] ; see also [10] . This method was used previously in the study of the diffusion coefficient 
for a Brownian particle in a periodic potential in |21j . Details about the numerical method can be 
found in Appendix [A} 

In all the numerical experiments we use a cosine potential, V{q) = Vq cos(uiiq), with u\ = 2ir/L. 
As a first test for the validity of our numerical method, in Figures [TJ [2] and [3] we compare the results 
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Figure 2: Solid Lines: D as a function of F for 7 = 0.01, 0.1 and 1. Markers: Monte Carlo 
estimates. Parameters of the simulations are as in Figure 1, with diffusivity now scaled by the value 
Dl = (/?7) _1 it would have in absence of the periodic potential 4>(x). 

obtained from the solution of the two PDEs with results obtained using Monte Carlo simulations. 
In particular, in Figures [I] and [2] we reproduce the results reported in [7] for 7 = 0.01 and go beyond 
this for larger values of 7. In all the Monte Carlo simulations reported in this paper we take a 
sufficiently large number of realizations, a sufficiently small time step and sufficiently long paths so 
that the results of the simulations are very accurate^] Details on the values of the parameters used 
in the simulations can be found in the caption figures. In Figure [3] we present results for U and D 
for larger values of 7. 

We emphasize the fact that the spectral method enables us to calculate the drift and diffusion 
coefficients very accurately for a very wide range of values of the friction coefficient 7 as well 
as the forcing F. As expected, the numerical method becomes computationally more expensive 
as 7 decreases, since more Hermite and Fourier modes are needed for the accurate calculation 
of the diffusion coefficient. We note also that, in two and higher dimensions, the underdamped 
regime requires appropriate preconditioning for the efficient solution of the resulting linear algebraic 
problem. 

1 In fact, in all the figures where the results of Monte Carlo simulations are presented, we also include the error 
bars. However, they are so small that they are barely visible. 
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a. U as a function of F for 7 = 0.5, 5 and 10. 



a 10 




b. D as a function of F for 7 = 0.5, 5 and 10. 



Figure 3: (a) Solid lines: U computed from (15). (b) Solid lines: D computed from (16a). 
Markers: Monte Carlo estimates. Parameters of the simulations are Vo = 1, (3 = 5, and L = 1. 
For the Monte Carlo simulations we have used a time step At = 0.01 over 6000000 time steps and 
averaging over 6000 trajectories. 
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10" 







0.2 



0.4 



0.6 



0.8 



Figure 4: U for 7 = 1 (blue lines) and 7 = 50 (red lines). Solid lines: U computed from 



(15). Markers: Monte Carlo estimates. Broken lines: Expansion for U in terms of F by solving 
numerically for the coefficients Vj given by (23a). 7 = 1: Dots: linear approximation, Dash: 
5th-order expansion, Dash-Dash: 9th-order expansion (overlapped with solid line). 7 = 50: 
Dots: linear approximation, Dash-dot: 3rd-order expansion, Dash-Dash: 5th-order expansion 
(overlapped with solid line). Parameters of the simulations are Vq = 1, /? = 5, and L = 1. For 
the Monte Carlo simulations we have integrated 5000 trajectories using a time step At = 0.01 over 
1000000 time steps for 7 = 1, while At = 0.005 over 4000000 time steps for 7 = 50. 
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Now we turn our attention to the numerical study of formulas (21 ) and (27). In Figure[4]we have 



calculated numerically the effective drift U using (15), and we have also calculated numerically the 



coefficients V n in (21). For this we need to solve the Poisson equations (24a), where the generator 
of the unperturbed dynamics, i.e. with F = 0, appears. We can see that as we increase the 
number of terms in the power series expansion, the series converges to the value of u computed 



from solving the stationary Fokker-Planck equation (14) and computing the integral in (15). We 



stress that, using the expansion (21) we can calculate the nonequilibrium drift for arbitrary values 



of the external forcing F using only information from the equilibrium dynamics. 

In Figure [5] we plot the diffusion coefficient D, as a function of the forcing F, using (16a). In 



addition, we plot the power series expansions of different orders according to assumption (27) that 



linear response relationships between drift and diffusion extend to the higher order coefficients. The 



drift U is computed as in Figure |4j using the expansion (21 ). While the power series expansion does 
match the value of D for F = 0, as it should according to linear response theory, the series does 
not converge to the values computed numerically for F 7^ using the spectral method described 
in Appendix [A} This shows in particular that the correction terms r.ke in Eq. (23c) are nontrivial, 
as also evidenced by the results of Figure 1 in [24J. To emphasize this point, we also compare 
in Figure [6] the effective diffusion coefficient with formula (27), where U is computed using (15) 
and the derivative with respect to F is approximated numerically using a centered finite difference 
scheme. See Figure [6j The linear response relations, however, do perform better for larger values 
of 7. 

Finally, we investigate the overdamped limit. The drift and diffusion coefficients of an over- 
damped particle moving in a one dimensional periodic potential under constant external force can 
be computed analytically in terms of quadratures. The exact formula for the effective drift is com- 
puted in ([29], [30, Ch. 9]), whereas the exact formula for the diffusion coefficient can be found 
in [21 and 1221. 



Expanding (14) and (17) in inverse powers of 7 we obtain 



and 



U 



D 



Uo 

7 



Do 

7 



+ 



+ 



<y<- 



(46) 



(47) 



where Uq and Dq denote the drift and diffusion coefficients for the overdamped problem (with 



7 scaled out, as described by the generator (49) below). Rather than computing the integrals 
in the formulas for Uq and Dq (see, e.g. [2^ 
Fokker-Planck and the Poisson equations |21j 



Eqn. B.6 and Eqn. B.9]) we solve the stationary 

Uo = f L (-V'(q) + F) po(q)dq, C oPo = 0, (48) 
Jo 

with £q the adjoint (Fokker-Planck operator) of the generator of the overdamped dynamics 

C = (-V'(q)+F)-? q+ r 1 -^- (49) 
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Figure 5: D for 7 = 1 (blue lines) and 7 = 50 (red lines). Solid lines :Homogenization for- 
mula (16a). Markers: Monte Carlo estimates. Broken lines: Expansion for D in terms of F 

1: Dots: 
50: 



assuming (27) and using the drift expansion (21). 7 



constant approximation, Dash: 
Dots: constant approximation, 



4th-order expansion, Dash-Dot: 8th-order expansion. 7 
Dash: 2nd-order expansion, Dash-Dot: 6th-order expansion. Parameters of the simulation are as 
in Figure |4} 
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Figure 6: D for 7=1 (blue lines) and 7 = 50 (red lines). Solid lines: Homogenization for- 
mula (16a). Broken lines: Computation of D assuming the Einstein relation (27), with derivative 
evaluated through centered finite differences of the drift computed from the homogenization for- 
mula (15). Parameters of the simulation as in Figure [4j 
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The generator is posed on [0, L] equipped with periodic boundary conditions. Similarly, the diffusion 
coefficient is given by 

Do = r 1 /(l + d q <(>o(q)) Po(q)dq, 

with 

-c (t>o(q) = -v\q) + F, 



on [0,L] with periodic boundary conditions. Higher order corrections in (46) and ( |47| ) can be 
obtained through the solution of further auxiliary Poisson equations. As shown in Figure [7] , the 
overdamped formulas for the drift and diffusion coefficients offer a very accurate approximation 
even for moderately high values of the friction coefficient, uniformly in F. 



5 Conclusions 

Using the framework of homogenization theory and multiscale analysis, we have developed a sys- 
tematic expansion of the effective drift and effective diffusivity for the nonequilibrium dynamics of 
a particle in a tilted periodic potential. The coefficients in this expansion relate the nonequilibrium 
transport coefficients to statistical averages involving the equilibrium dynamics (with no imposed 
tilt), computed through the solutions of boundary value problems for deterministic partial differ- 
ential equations of hypoelliptic type. The expansions give a detailed description of how Einstein's 
relation between the diffusivity and mobility of a particle is violated in higher orders with respect to 
the perturbation from equilibrium. Our theoretical results were confirmed by numerical simulations 
based on a new efficient spectral method for the solution of Poisson equations for the generator of 
the Langevin dynamics. 

Our method of analysis can be readily extended, with suitable elaboration of notation, to mul- 
tiple dimensions. Other substantial directions for future research include the application of the ho- 
mogenization procedure to multiscale and locally periodic potentials, as well as to time-dependent 
external forcing. This last setting could have particular relevance to the study of stochastic reso- 
nance phenomena. 
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A Numerical Algorithm 

In this appendix we present a numerical approach for solving the 1-dimensional stationary Fokker- 
Planck equation (14) together with the cell problem (17) for computing V and D via (15) and (16a). 
This numerical method is based on the approach by [27J and consists in a spectral decomposition 
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Figure 7: (a) Solid lines: Effective drift U computed from (15). Broken lines: Asymptotic 



approximation (46). (b) Solid lines: Effective diffusivity D computed from (16a). Broken lines: 



Asymptotic approximation (47). Parameters of the simulations are Vq = 1, /3 = 5, and L = 1. 
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of the solution of (14) and (17) in terms of Hermite polynomials and Fourier series, followed by 



a recursive method to solve the resulting system of algebraic equations. Since this approach is 
presented in pTj for finding numerically pp{q,p) and U , we will focus on the computation of D via 



the solution for the auxiliary field 4>(q,p) in equation ((|17j)) and equation (16a). 



A.l Solution in terms of Hermite polynomials. 

The cell problem for the auxiliary field (f>(q,p) can be written in terms of the infinitesimal generator 
of the Ornstein-Uhlenbeck (OU) process as introduced in Section [3j 

S = -pd p + f3- 1 d 2 p , 

as 

£4>(q,p) = pd q( f> + (-V'(q) + F)d p ct> + 1 S^=U-p, 



with U the effective drift as given by (15). Note that the invariant distribution p(p) of the OU 



process, S* p(p) = 0, implies p ~ e I 2 . In view of the structure of the previous equation, we use 
the following representation for its solution, 



4>{q, p) = ^2 < t ) n{q)H n {p), 



(50) 



n=0 



where 4> n {q) is a series of functions to be determined. H n (p) are rescaled Hermite polynomials 



H n (p) 

Be n (x) 



He n (p/3 



1/2N 



n r x 2 /2 d n c -x 2 /2 

dx 11 



(-l) n e~ e 



which are the eigenfunctions of the operator S, 

SH n (p) = -nH n (p), n = l,2, ... 

Also, these rescaled Hermite polynomials are orthonormal with respect to the unperturbed station- 
ary distribution: 

/oo 
H n (p)H m (p)p(p) dy = 5 nm , 
-oo 

and satisfy the relations : 

P H n (p) = /T 1 / 2 (Vn+TH re+ i(p) + V^H n -i(p)) , 
H' n (p) = {Pn) l/2 H n -i{p). 



Upon substituting (50) into ((17)), projecting against Ho, Hi, and H n for n > 2 respectively 



and using the orthonormality property of the Hermite polynomials, we obtain the following infinite 
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system of ordinary differential equations for (f> n (q), 

( f ) [( q ) + f](-V / (q) + F)MQ) = VPU, 
0d(g) " iVPMv) + P {-V'(q) + F) V2MQ) + V2<j>' 2 (q) = -1, 
V"Ci(9) ~ l\fj3n<t>n(q) + VnTl 4>' n+ i{q) 

+P (-V'(q) + F) V^+l <An+l (q) = 0, 

for n = 2, 3, . . . 

A. 2 Spectral decomposition. 

Since the solution to the cell problem must be periodic in q, the auxiliary functions <f> n (q) must 
also be periodic. It is natural then to express these functions in terms of their Fourier series, 



2vrj 

j=-oo 



L 

For simplicity, we will focus now on the simplest periodic potential, namely, 

V{q) = F cos(2^q/L), 

although more complex potentials can be studied. In terms of this potential, the equations take 
the following form, 

iuM+P - + F*{\ = y/0U6 jja , 

iu) j & -j<yi3<f>{+ 

V^(^< + i + /3(^(C\-<tl)+^l + i)) = 0, (51) 

for n = 2,3,..., j = 0, ±1, ±2, . . . 

A. 3 Solution of <p n . 

We now proceed to describe the numerical algorithm for computing D. In order to solve (17) in its 



spectral representation (51), we approximate 4> n (q) by a Galerkin truncation of the Fourier series 
after the Mth term, 

M 

Ml)* E ®n^ q - 
j=-M 



The infinite system of algebraic equations (51) becomes then an infinite, tri-diagonal system of 



equations expressed as follows. By explicitly writing the real and imaginary parts of = £n + iff 
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and using the fact that the solution must be real- valued (which implies that ^ n 3 = £n and r} n 3 = —rjn 



) we form the vectors, 



6 J1 

Sri 



ill 



, 71 = 0,1,2, 



This representation leads to the following system of equations, 

Q+* n _i + Q n *„ + Q-* n+ i = 0, n = 2,3,... 
These matrices are given, for n = 0, 1, . . ., by, 

Qn = -1 ^fPnI 2 M+l, 

where is the k x k identity matrix. For n = 1, 2, ... we have, 

/ | 

| -wi 





... 



wi 
uj 2 

^ 



... 
... 

UJ M 



o \ 




(52a) 
(52b) 
(52c) 



Qaa Qab 



Q 



ba 



Q 



bb 



Qaa = F f3 Im+1, Qbb = F P Im 
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Q 



ab 



( -PV uji 

-wi -pVoUi/2 

PV u)x/2 -oj 2 -PVoUi/2 









V 







/3F wi/2 -u M J 



( -pV u>i/2 




Q 



bit 



ui 

-PVoUi/2 u 2 PVoUi/2 



\ 




V 



-pV ui/2 u M ) 



[B] k = y/pU6 ktli [A] k = -5 k>1 , 

where [B] k represents the fcth element of the vector B (respectively for A.) In order to solve the infi- 
nite system of algebraic equations, we impose some boundary condition of the form <&at+i = Sn*&n, 
for large N. Tested boundary conditions include Sn = (Dirichlet boundary condition), which 
we employed in the simulations in Section [4] , and Sn = I2M+1 (Neumann boundary condition). 
Defining matrices {S n }n=Q recursively downwards from n = N — 1 by 

S n = -(Q n+l + Q- +1 S n+1 )~ 1 Q+ +1 for n = 0, . . . , N - 1, 

we can check by induction (again downwards) that for n = 1, . . . , N, 

$ n+ i = S n & n (53) 

Indeed, this relation is already in force for n = N, and assuming it to be true for some n = m > 2, 



from Eq. (52c) we find: 

Q^m-l + Qm®m + Q~S m ^m = Qm*m-1 + [Q m + Q m S m ) *m = 0, 



so that Eq. ( 53 ) holds for n = m — 1 as well. Turning now to Eqs. fl52aj ) and ( |52b[ ), we have 

Q+*o + (Qi + Qr^i)*i = A (54a) 
Q *i = S > ( 54b ) 
from which we find by solving for 3>i in terms of 3>o : 

^ 1 = S ^ Q +{Q 1 + QjSi)~ 1 A. 



Substituting this expression into Eq. (54b), we finally obtain a closed equation for 3>o : 

-i 



Q S 3>o = B-Q (Qi + QrSi) A. 



(55) 
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The matrix Q Q So will have one null eigenvalue (corresponding to the null space of £). One can 
verify, by considering the analogous numerical solution scheme for pa and U, presented in |27j . 



that the right hand side of Eq. (55) satisfies the solvability condition that it be orthogonal to 



the left eigenvector of Q Sq with zero eigenvalue. A unique solution for 3>o is then obtained by 



discretization of the auxiliary condition in Eq. (17). In particular, representing the solution to the 



stationary Fokker-Planck ( 14 ) by a Hermite polynomial expansion 



pp(q,p) = p(p) Rn(q)H n (p), 



n=0 



and approximating the functions R n (q) by a finite Fourier series, with coefficients organized into 



vectors R n analogously to Eq. (50), this auxiliary condition reads: 



N 
n=0 



This then determines, with Eq. (55), <&o from which he remaining {*& n }n=i are found recursively 



using the matrices S n and the relations (53). Once the vectors <&j are found, D is easily computed 



by replacing the proposed solution for pp{q,p) and (f)(q,p) in ( 16a) and using the Hermite polynomial 
properties to obtain: 



N 



D 



n=0 



n + 1 



[R T n+l ®n]i+2R T n 3> n+1 



n+l]lj 



B Alternative Approach to Obtaining Corrections to Einstein's 
Formula 



The relation between the diffusivity and mobility is expressed in [2] as follows (in our notation): 

(56) 



D = P~ ld ^+ lim -i 

dF T-s>oo 27 j 



T 



where (g; h) = {gh) — (g){h) and (•) denotes an average over the stochastic noise (and possibly 
random initial conditions). The correction term was studied in [2] on the model system as well 
as other non-equilibrium systems through direct numerical simulation of the governing dynamical 



equations and Monte Carlo estimation of the statistical average. We can express Eq. (56) in terms 



of deterministic operators through the following formal manipulations. First, we re-express 

q(T)-q(0)= f T p(t')dt', 
Jo 

which avoids the complication of working with the nonperiodic variable q(t). We then have: 

D = r 1 ^+lim— / / (p(t');-V'( q (t)) + F)dt'dt 
dF T^oo 2jT Jo Jo ' 



25 



Now, thanks to the large factor of T in the denominator, we may neglect initial transients and 
evaluate the statistical average in the nonequilibrium steady state, i.e., with single-time statistics 
governed by the invariant density pp, the solution of the stationary Fokker-Planck equation (14). 
We then express the two-time correlation function formally using the evolution operator e - ^*, 
where C denotes the generator of the Langevin dynamics, for the forward-in-time variable, and the 
projection operator 

Pg = g- (g) P 



where 



to obtain: 



/ gppdpdq 
Jx 



D 



P 



^dU 
dF 



lim — - , 

T^oo 27T Jo 

T 

e 



(Pp) (e c ^PV'(q)^ dt' 



+ 



C(f-t) ¥p \ (f> V '(q))) dt' 



dt 



dU 



lim 



1 



£- 2 {e CT - I - CT)Pp) (PV'{q))) 



dF T^oo 2jT 

+ ((Pp) (£- 2 (e CT - I - CT)FV'(q))] 



where I is the identity operator. Using now the nonpositivity of the operator C and the fact that 
RanP = Ran£ since P projects onto the subspace orthogonal to the kernel of £*, the 1? adjoint 
of C, we can evaluate the T — > 00 limit to obtain the following formal operator-theoretic equivalent 
to the formula Eq. ( 56 ) from [2] : 



D = p 



-xdU 



l r 



dF + 2 7 



[{C-'Fp) V'(q)) + (QPp) (jC-Wfa))) 



(57) 



A somewhat more compact formula can be obtained by defining the L 2 (p) adjoint of £, which can 
be computed as: 

C = C -Fa- +2j{p + p- 1 d p ln Pl3 )a-, (58) 
where Cq is defined at the end of Proposition 3.1, so that 



D 



(59) 



In both of these expressions, we note that 



p — U, from Eq. ( 15 ) . 



Inspecting expression Eq. (57) for the correction to the Einstein relation, we see that beyond 
computing pp as the stationary solution of the Fokker-Planck equation (14), we must solve a Poisson 
equation of the form (17) as well as a second Poisson equation of the form 

-Ciji = PV'. 



In the expression Eq. (59), we must solve a stationary Fokker-Planck equation (14), a Poisson 
equation of the form (17), as well as an adjoint-Poisson equation of the form 

-Li) = Pp. 
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In both cases, it seems that an additional equation would need to be solved beyond the stationary 
Fokker-Planck equation (14) and a single Poisson equation (17) necessary in the homogenization 
approach. On the other hand, computing the mobility at general values of the tilt F from the 
nonperturbative homogenization equations would require a differentiation between different values 
of F. The direct formula (56) would generally of course need to be evaluated through Monte Carlo 
averages involving a large number of sample trajectories run for sufficiently long. 

The perturbation theory with respect to F developed for the homogenization equations in 
Section [3] has the virtue of allowing the simultaneous numerical computation of the diffusivity and 
drift for a range of values of tilt F, rather than one value at a time. One could introduce similar 
perturbation expansions with respect to tilt F into the formulas (57) and (59). We attempted 
to examine whether this would give equivalent results, but found this effort frustrating. On the 
one hand, computing Eq. (57) perturbatively would introduce the perturbative series solution to 
a second Poisson equation completely absent from the homogenization theory, so it would be 
difficult to relate the results. Expression (59) has more promise because to leading order, 
is identical to the simple operator Cq , which is just a time reversal of the operator Cq 1 . However, 



implementing the perturbation expansion on Eq. (59), even to first order, produced considerably 
more unwieldy equations than emerged from the homogenization equations, and again how to 
relate the resulting expressions was unclear. The main complication is the propagation of the 
perturbation expansion (32) for the invariant density through the adjoint operator C (58). Perhaps 
a more clever analysis would provide a linkage between the formula for the correction (56) to the 
Einstein relation from [2] and the perturbative expansion we have developed in Proposition 3.1 but 
it appears that computations are considerably simpler by conducting the perturbation expansion 
on the homogenization equations as we have done in Section [3j 
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